I fixed the TPCF issues, but now the delta sigma projection integration is being troublesome. I'm gonna do them for a few rp_bins to figure out what is sensible agian.


In [1]:
import numpy as np
from glob import glob
from os import path

In [2]:
from matplotlib import pyplot as plt
%matplotlib inline
import seaborn as sns
sns.set()
import matplotlib.colors as colors

In [3]:
from pearce.mocks.kittens import DarkSky
from pearce.mocks import tpcf as pearce_tpcf
from halotools.mock_observables import tpcf as halotools_tpcf
from halotools.empirical_models import Zheng07Cens, Zheng07Sats
from collections import OrderedDict
from time import time
from scipy.optimize import minimize_scalar
import yaml

In [4]:
output_dir = './'# '/home/users/swmclau2/Git/pearce/bin/covmat/ds14_covmat/'

In [5]:
config_fname = 'xi_cosmo_trainer.yaml'

with open(path.join(output_dir, config_fname), 'r') as ymlfile:
    cfg = yaml.load(ymlfile)

nd = float(cfg['HOD']['fixed_nd'] )
min_ptcl = int(cfg['HOD']['min_ptcl'])
r_bins = np.array(cfg['observation']['bins'] ).astype(float)

hod_param_ranges =  cfg['HOD']['ordered_params']


logMmin_bounds = hod_param_ranges['logMmin']


del hod_param_ranges['logMmin']

In [6]:
def make_LHC(ordered_params, N, seed = None):

    if seed is None:
        seed = int(time())
    np.random.seed(seed)

    points = []
    # by linspacing each parameter and shuffling, I ensure there is only one point in each row, in each dimension.
    for plow, phigh in ordered_params.itervalues():
        point = np.linspace(plow, phigh, num=N)
        np.random.shuffle(point)  # makes the cube random.
        points.append(point)
    return np.stack(points).T


def add_logMmin(hod_params, cat):

    hod_params['logMmin'] = 13.0 #initial guess
    #cat.populate(hod_params) #may be overkill, but will ensure params are written everywhere
    def func(logMmin, hod_params):
        hod_params.update({'logMmin':logMmin})
        return (cat.calc_analytic_nd(hod_params, min_ptcl = min_ptcl) - nd)**2

    res = minimize_scalar(func, bounds = logMmin_bounds, args = (hod_params,), options = {'maxiter':100}, method = 'Bounded')

    # assuming this doens't fail
    #print 'logMmin', res.x
    hod_params['logMmin'] = res.x

In [7]:
from pearce.mocks.kittens import TestBox
cat = TestBox(boxno = 0, realization = 0, system = 'sherlock')
cat.load(1.0, HOD='zheng07', particles = True, downsample_factor = 1e-2)

In [8]:
# TODO seed here for constant HODs
# TODO maybe just do 5, 10 may be overkill
N = 10
LHC = make_LHC(hod_param_ranges, N, 24)
hod_dicts = [dict(zip(hod_param_ranges.keys(), vals)) for vals in LHC]

In [9]:
cat.populate(hod_dicts[1])

In [10]:
cat._get_cosmo_param_names_vals()


Out[10]:
(['ombh2', 'omch2', 'w0', 'ns', 'ln10As', 'H0', 'Neff'],
 array([  2.32629000e-02,   1.07830000e-01,  -7.26513000e-01,
          9.80515000e-01,   3.03895000e+00,   6.32317000e+01,
          2.95000000e+00]))

In [11]:
from halotools.mock_observables import return_xyz_formatted_array, delta_sigma, tpcf

In [12]:
def calc_ds(cat, bins):
    #n_cores = self._check_cores(n_cores)
    n_cores = 2
    
    x_g, y_g, z_g = [cat.model.mock.galaxy_table[c] for c in ['x', 'y', 'z']]
    pos_g = return_xyz_formatted_array(x_g, y_g, z_g, period=cat.Lbox)

    x_m, y_m, z_m = [cat.halocat.ptcl_table[c] for c in ['x', 'y', 'z']]
    pos_m = return_xyz_formatted_array(x_m, y_m, z_m, period=cat.Lbox)

    rp_bins = bins #if not angular else self._rp_from_ang(bins)

    # Halotools wnats downsampling factor defined oppositley
    # TODO verify little h!
    # TODO maybe split into a few lines for clarity
    return delta_sigma(pos_g / cat.h, pos_m / cat.h, cat.pmass / cat.h,
                   downsampling_factor=1. / cat._downsample_factor, rp_bins=rp_bins,
                   period=cat.Lbox / cat.h, num_threads=n_cores, cosmology=cat.cosmology)[1] / (1e12)

In [13]:
from scipy.interpolate import interp1d
#import pyccl as ccl
from scipy.integrate import quad
def calc_ds_analytic(cat, bins, xi = None, rbins = None): n_cores = 1#self._check_cores(n_cores) # calculate xi_gg first xi_kwargs = {} if xi is None: assert rbins is None _rbins = np.logspace(-2, 1.6, 21) rbins = np.zeros((_rbins.shape[0]-1)) rbins[0] = _rbins[0] rbins[1:] = _rbins[2:] xi = cat.calc_xi_gm(rbins, n_cores=n_cores, **xi_kwargs) else: assert rbins is not None, "Must specify rbins along with xi" if np.any(xi <= 0): warnings.warn( "Some values of xi are less than 0. Setting to a small nonzero value. This may have unexpected behavior, check your HOD") xi[xi <= 0] = 1e-3 rpoints = (rbins[:-1] + rbins[1:]) / 2.0 xi_rmin, xi_rmax = rpoints[0], rpoints[-1] # make an interpolator for integrating xi_interp = interp1d(np.log10(rpoints), np.log10(xi)) # get the theotertical matter xi, for large scale estimates names, vals = cat._get_cosmo_param_names_vals() param_dict = {n: v for n, v in zip(names, vals)} if 'omch2' in param_dict: # in other units, convert new_param_dict = {} new_param_dict['Omega_c'] = param_dict['omch2']/cat.h**2 new_param_dict['Omega_b'] = param_dict['ombh2']/cat.h**2 new_param_dict['n_s'] = param_dict['ns'] new_param_dict['h'] = cat.h new_param_dict['A_s'] = np.exp(param_dict['ln10As'])/(np.power(10, 10)) param_dict = new_param_dict elif 'Omega_c' not in param_dict: param_dict['Omega_c'] = param_dict['Omega_m'] - param_dict['Omega_b'] del param_dict['Omega_m'] cosmo = ccl.Cosmology(**param_dict) big_rbins = np.logspace(1, 2.3, 21) big_rpoints = (big_rbins[1:] + big_rbins[:-1]) / 2.0 big_xi_rmax = big_rpoints[-1] xi_mm = ccl.correlation_3d(cosmo, cat.a, big_rpoints) xi_mm[xi_mm < 0] = 1e-6 # may wanna change this? xi_mm_interp = interp1d(np.log10(big_rpoints), np.log10(xi_mm)) # correction factor bias = np.power(10, xi_interp(1.2) - xi_mm_interp(1.2)) rhocrit = cat.cosmology.critical_density(0).to('Msun/(Mpc^3)').value rhom = cat.cosmology.Om(0) * rhocrit * 1e-12 # SM h^2/pc^2/Mpc; integral is over Mpc/h def sigma_integrand_medium_scales(lRz, Rp, xi_interp): Rz = np.exp(lRz) #print "Med", Rz, Rp, rpoints[0] return Rz * 10 ** xi_interp(np.log10(Rz * Rz + Rp * Rp) * 0.5) def sigma_integrand_large_scales(lRz, Rp, bias, xi_mm_interp): Rz = np.exp(lRz) #print "Large", Rz, Rp, rpoints[0] return Rz * bias * 10 ** xi_mm_interp(np.log10(Rz * Rz + Rp * Rp) * 0.5) ### calculate sigma first### sigma_rpoints = np.logspace(-1.5, 2.2, 15) sigma = np.zeros_like(sigma_rpoints) for i, rp in enumerate(sigma_rpoints): log_u_ss_max = np.log(xi_rmax ** 2 - rp ** 2) / 2.0 # Max distance to integrate to log_u_ls_max = np.log(big_xi_rmax ** 2 - rp ** 2) / 2.0 # Max distance to integrate to if not np.isnan(log_u_ss_max): # rp > xi_rmax small_scales_contribution = \ quad(sigma_integrand_medium_scales, -10, log_u_ss_max, args=(rp, xi_interp))[0] Rz = np.exp(log_u_ls_max) large_scales_contribution = quad(sigma_integrand_large_scales, log_u_ss_max, log_u_ls_max, \ args=(rp, bias, xi_mm_interp))[0] elif not np.isnan(log_u_ls_max): small_scales_contribution = 0 large_scales_contribution = quad(sigma_integrand_large_scales, np.log(xi_rmax), log_u_ls_max, \ args=(rp, bias, xi_mm_interp))[0] else: small_scales_contribution = large_scales_contribution = 0.0 assert not any(np.isnan(c) for c in (small_scales_contribution, large_scales_contribution)), "NaN found, aborting calculation" sigma[i] = (small_scales_contribution + large_scales_contribution) * rhom * 2; sigma_interp = interp1d(np.log10(sigma_rpoints), sigma) ### calculate delta sigma ### def DS_integrand_medium_scales(lR, sigma_interp): #print 'DS', np.exp(lR), rp_points[0] return np.exp(2 * lR) * sigma_interp(np.log10(np.exp(lR))) rp_bins = bins #if not angular else cat._rp_from_ang(bins) rp_points = (rp_bins[1:] + rp_bins[:-1]) / 2.0 lrmin = np.log(sigma_rpoints[0]) ds = np.zeros_like(rp_points) for i, rp in enumerate(rp_points): result = quad(DS_integrand_medium_scales, lrmin, np.log(rp), args=(sigma_interp,))[0] print i, np.exp(lrmin), rp, result * 2 / (rp ** 2), sigma_interp(np.log10(rp)) # print result, rp, sigma_interp(np.log10(rp)) ds[i] = result * 2 / (rp ** 2) - sigma_interp(np.log10(rp)) return xi, sigma, ds

In [14]:
from halotools.mock_observables import delta_sigma_from_precomputed_pairs
from halotools.mock_observables import  total_mass_enclosed_per_cylinder as tmepc_halotools
from halotools.mock_observables.surface_density.mass_in_cylinders import _enclosed_mass_process_args
from halotools.mock_observables.surface_density.weighted_npairs_per_object_xy import weighted_npairs_per_object_xy

In [15]:
def total_mass_enclosed_per_cylinder(centers, particles,
        particle_masses, downsampling_factor, rp_bins, period,
        num_threads=1, approx_cell1_size=None, approx_cell2_size=None):

#  Perform bounds-checking and error-handling in private helper functions
    #print period
    args = (centers, particles, particle_masses, downsampling_factor,
        rp_bins, period, num_threads)
    result = _enclosed_mass_process_args(*args)
    centers, particles, particle_masses, downsampling_factor, \
        rp_bins, period, num_threads, PBCs = result

    #print 'PBCs:', PBCs
    #print 'Period:', period
    #period = None
    mean_particle_mass = np.mean(particle_masses)
    normalized_particle_masses = particle_masses/mean_particle_mass

    # Calculate M_tot(< Rp) normalized with internal code units
    total_mass_per_cylinder = weighted_npairs_per_object_xy(centers, particles,
        normalized_particle_masses, rp_bins,
        period=None, num_threads=num_threads, #try large finite PBCs
        approx_cell1_size=approx_cell1_size,
        approx_cell2_size=approx_cell2_size)

    # Renormalize the particle masses and account for downsampling
    total_mass_per_cylinder *= downsampling_factor*mean_particle_mass

    return total_mass_per_cylinder

In [16]:
def calc_tm(cat, bins):
    #n_cores = self._check_cores(n_cores)
    n_cores = 4
    
    x_g, y_g, z_g = [cat.model.mock.galaxy_table[c] for c in ['x', 'y', 'z']]
    pos_g = return_xyz_formatted_array(x_g, y_g, z_g, period=cat.Lbox)

    x_m, y_m, z_m = [cat.halocat.ptcl_table[c] for c in ['x', 'y', 'z']]
    pos_m = return_xyz_formatted_array(x_m, y_m, z_m, period=cat.Lbox)

    rp_bins = bins #if not angular else self._rp_from_ang(bins)
    tm = total_mass_enclosed_per_cylinder(pos_g/cat.h, pos_m/cat.h, cat.pmass/cat.h,\
                                     1./cat._downsample_factor, rp_bins,# cat.Lbox/cat.h,\
                                          cat.Lbox/cat.h,\
                                         num_threads = n_cores)
    print 
    tm_ht = tmepc_halotools(pos_g/cat.h, pos_m/cat.h, cat.pmass/cat.h,\
                                     1./cat._downsample_factor, rp_bins,# cat.Lbox/cat.h,\
                                          cat.Lbox/cat.h,\
                                         num_threads = n_cores)
    
    return tm, tm_ht

In [17]:
def calc_ds2(cat, bins, tm):
    
    x_g, y_g, z_g = [cat.model.mock.galaxy_table[c] for c in ['x', 'y', 'z']]
    pos_g = return_xyz_formatted_array(x_g, y_g, z_g, period=cat.Lbox)

   
    ds = delta_sigma_from_precomputed_pairs(pos_g / cat.h, tm, bins, \
                                              cat.Lbox/cat.h, cosmology=cat.cosmology)[1] / (1e12)
    
    return ds

In [18]:
rp_bins = np.logspace(-0.5, 2.0, 11)

In [19]:
names, vals = cat._get_cosmo_param_names_vals()
param_dict = {n: v for n, v in zip(names, vals)}

new_param_dict = {}
new_param_dict['Omega_c'] = param_dict['omch2']/cat.h**2
new_param_dict['Omega_b'] = param_dict['ombh2']/cat.h**2
new_param_dict['n_s'] = param_dict['ns']
new_param_dict['h'] = cat.h
new_param_dict['A_s'] = np.exp(param_dict['ln10As'])/(np.power(10, 10))

In [20]:
param_dict


Out[20]:
{'H0': 63.231699999999996,
 'Neff': 2.9500000000000002,
 'ln10As': 3.0389499999999998,
 'ns': 0.98051499999999991,
 'ombh2': 0.023262900000000003,
 'omch2': 0.10783,
 'w0': -0.72651299999999996}

In [21]:
new_param_dict


Out[21]:
{'A_s': 2.0883304249678889e-09,
 'Omega_b': 0.058182735712595794,
 'Omega_c': 0.26969313335350292,
 'h': 0.632317,
 'n_s': 0.98051499999999991}

In [22]:
rp_bins


Out[22]:
array([   0.31622777,    0.56234133,    1.        ,    1.77827941,
          3.16227766,    5.62341325,   10.        ,   17.7827941 ,
         31.6227766 ,   56.23413252,  100.        ])
xi, sigma, ds_analytic = calc_ds_analytic(cat, rp_bins)
ds = calc_ds(cat, rp_bins)
print ds

In [23]:
tm, tm_ht = calc_tm(cat, rp_bins)




In [24]:
ds = calc_ds2(cat, rp_bins, tm)

In [26]:
print ds


[ 22.67152071  12.91536306   6.73179328   3.23496689   1.59690043
   0.94686751   0.71960138   0.6919193    0.83918369   1.24496502]

In [28]:
rp_bins


Out[28]:
array([  0.31622777,   0.56234133,   1.        ,   1.77827941,
         3.16227766,   5.62341325,  10.        ,  17.7827941 ,
        31.6227766 ,  56.23413252, 100.        ])
for idx in xrange(10): plt.hist(np.log10(tm[:, idx] + 1e-6), bins = np.linspace(12, 18, 100) , alpha = 0.5); plt.hist(np.log10(tm_ht[:,idx] + 1e-6) , bins = np.linspace(12, 18, 100), alpha = 0.5); #plt.loglog() plt.show();

In [42]:
tm.shape


Out[42]:
(796509, 11)

In [77]:
idxs = cat.model.mock.galaxy_table['x'] < 10#, \
                    # cat.model.mock.galaxy_table['x'] > 999.99)

In [80]:
plt.hist(np.log10(tm[idxs,idx] + 1e-6) , bins = 100, alpha = 0.3, label = 'No PBC');
plt.hist(np.log10(tm_ht[idxs,idx] + 1e-6) , bins = 100, alpha = 0.3, label = 'PBC');



In [79]:
plt.hist(np.log10(tm[~idxs,idx] + 1e-6) , bins = 100, alpha = 0.3);
plt.hist(np.log10(tm_ht[~idxs,idx] + 1e-6) , bins = 100, alpha = 0.3);



In [49]:
plt.hist(np.log10(tm[:,idx] + 1e-6) , bins = 100, alpha = 0.3);
plt.hist(np.log10(tm_ht[:,idx] + 1e-6) , bins = 100, alpha = 0.3);


Out[49]:
(array([3.0000e+00, 5.0000e+00, 1.0000e+00, 6.0000e+00, 7.0000e+00,
        5.0000e+00, 5.0000e+00, 7.0000e+00, 4.0000e+00, 9.0000e+00,
        1.0000e+01, 9.0000e+00, 8.0000e+00, 1.0000e+01, 1.2000e+01,
        1.6000e+01, 1.6000e+01, 1.0000e+01, 1.7000e+01, 2.2000e+01,
        1.7000e+01, 2.3000e+01, 2.1000e+01, 1.9000e+01, 2.6000e+01,
        3.1000e+01, 2.9000e+01, 3.5000e+01, 3.2000e+01, 5.1000e+01,
        4.8000e+01, 5.2000e+01, 7.3000e+01, 8.5000e+01, 1.1600e+02,
        1.7200e+02, 2.5400e+02, 3.2700e+02, 4.7000e+02, 5.9600e+02,
        7.1400e+02, 8.3700e+02, 9.4100e+02, 1.0300e+03, 1.0870e+03,
        1.1870e+03, 1.2830e+03, 1.2760e+03, 1.4350e+03, 1.3670e+03,
        1.3790e+03, 1.3930e+03, 1.4710e+03, 1.4870e+03, 1.5190e+03,
        1.5010e+03, 1.5990e+03, 1.5610e+03, 1.6280e+03, 1.6540e+03,
        1.6980e+03, 1.6800e+03, 1.7630e+03, 1.7520e+03, 1.7940e+03,
        1.9240e+03, 1.8970e+03, 1.9280e+03, 1.9770e+03, 2.0220e+03,
        2.0450e+03, 2.1780e+03, 2.2840e+03, 2.3380e+03, 2.4530e+03,
        2.4300e+03, 2.5000e+03, 2.9410e+03, 4.0170e+03, 6.0910e+03,
        1.3087e+04, 2.5209e+04, 3.7993e+04, 4.6703e+04, 6.1439e+04,
        7.9624e+04, 8.4397e+04, 9.1546e+04, 7.5861e+04, 6.9424e+04,
        5.2046e+04, 3.6217e+04, 2.1991e+04, 1.1594e+04, 6.5350e+03,
        2.8270e+03, 1.0470e+03, 1.0060e+03, 1.0850e+03, 1.5800e+02]),
 array([17.21073651, 17.21726689, 17.22379726, 17.23032764, 17.23685801,
        17.24338839, 17.24991876, 17.25644914, 17.26297951, 17.26950989,
        17.27604026, 17.28257064, 17.28910102, 17.29563139, 17.30216177,
        17.30869214, 17.31522252, 17.32175289, 17.32828327, 17.33481364,
        17.34134402, 17.34787439, 17.35440477, 17.36093515, 17.36746552,
        17.3739959 , 17.38052627, 17.38705665, 17.39358702, 17.4001174 ,
        17.40664777, 17.41317815, 17.41970852, 17.4262389 , 17.43276928,
        17.43929965, 17.44583003, 17.4523604 , 17.45889078, 17.46542115,
        17.47195153, 17.4784819 , 17.48501228, 17.49154265, 17.49807303,
        17.5046034 , 17.51113378, 17.51766416, 17.52419453, 17.53072491,
        17.53725528, 17.54378566, 17.55031603, 17.55684641, 17.56337678,
        17.56990716, 17.57643753, 17.58296791, 17.58949829, 17.59602866,
        17.60255904, 17.60908941, 17.61561979, 17.62215016, 17.62868054,
        17.63521091, 17.64174129, 17.64827166, 17.65480204, 17.66133242,
        17.66786279, 17.67439317, 17.68092354, 17.68745392, 17.69398429,
        17.70051467, 17.70704504, 17.71357542, 17.72010579, 17.72663617,
        17.73316655, 17.73969692, 17.7462273 , 17.75275767, 17.75928805,
        17.76581842, 17.7723488 , 17.77887917, 17.78540955, 17.79193992,
        17.7984703 , 17.80500067, 17.81153105, 17.81806143, 17.8245918 ,
        17.83112218, 17.83765255, 17.84418293, 17.8507133 , 17.85724368,
        17.86377405]),
 <a list of 100 Patch objects>)

In [27]:
rpoints = (rp_bins[1:] + rp_bins[:-1])/2.0
plt.plot(rpoints, ds.squeeze(), label = 'Halotools')
#plt.plot(rpoints, ds2.squeeze(), label = 'Halotools 2')
#plt.plot(rpoints, ds2_ht.squeeze(), label = 'Halotools 3')

#plt.plot(rpoints, ds_analytic.squeeze(), label = 'Analytic')

plt.legend(loc='best')
plt.loglog();



In [ ]:
ds/ds2

In [ ]:
ds/ds2_ht

To test: set period to large finite value set all periods to none run with larger rp bins to see if problems continue at larger scaels